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The DAMA experiments are running deep underground in the Gran Sasso National Laboratory. Several 
interesting results have been achieved so far. Here a maximum likelihood method to search for the WIMP 
annual modulation signature is discussed and applied to a set of preliminary test data collected with large mass 
highly radiopure Nal(Tl) detectors. Various related technical arguments are briefly addressed. 



1. Introduction. 

The DAMA experiments have already achieved 
several interesting results using various target- 
detectors |-|]. This paper is devoted to discuss 
an analysis strategy on annual modulation signa- 
ture; a practical example is performed by using 
a preliminary data set collected with large-mass 
highly radiopure Nal(Tl) detectors. Standard as- 
sumptions for the WIMP model have been con- 
sidered. A more extensive discussion on analysis 
strategies is in 

In the direct searches by WIMP-nucleus elas- 
tic scattering, the possible presence of a WIMP 
signal can be extracted from the background by 
the annual modulation of the WIMP rate in the 
target-detector ^,^-0]. In fact, the expected 
recoil energy spectrum depends on the WIMP 
velocity distribution and on the Earth velocity 
in the galactic frame, v r (t), which varies along 
the year according to the expression: v r (t) = 

Vs un + V E arthCOS"fCOSUj(t - t ) . Here V S un = 232 

km/s is the Sun velocity with respect to the galac- 
tic halo and V Earth — 30 km/s is the Earth orbital 
velocity around the Sun on a plane with inclina- 



tion 7 = 60° with respect to the galactic one; fur- 
thermore, ui — 27r/T with T=l year and to ~ 2 nd 
June (when the Earth speed is at maximum). The 
WIMP velocity distribution in the galactic halo 
frame is considered to be a Maxwellian distribu- 
tion with vo parameter (defined as \/2/3 • v r , m . s ,) 
and a cut-off velocity equal to the escape veloc- 
ity. The Earth velocity can be conveniently ex- 
pressed in unit of vo : 77(f) = v r (t)/v = r/a + 
Arj cosuj(t — to), being ij =1.05 the yearly aver- 
age of n and An = 0.07. Since Arj « 770, the 
signal rate in the k-th energy interval is accu- 
rately given by the first order Taylor approxima- 
tion: S k [v(t)} = S k [r/ } + [^} Vo AriCosu(t-t ) = 
So,k+S mt kCosuj(t—to). The contribution from the 
highest order terms is lower than 0.1%. To select 
the S m .k the highest sensitivity is obtained when 
considering the smallest energy bins allowed by 
the available statistics and whole year data tak- 
ing Q). Although the time dependent effect is 
expected to be small, a suitable large-mass low- 
radioactive set-up with an efficient stability mon- 
itoring can point out its presence by time corre- 
lation analysis Note that, with the present 
technology, the annual modulation remains the 
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main signature of a possible WIMP signal]]. 

2. The experimental data 

The data considered here have been collected 
with nine 9.70 kg Nal(Tl) detectors, part of the 
115.5 kg highly radiopure Nal(Tl) set-up now 
running at the Gran Sasso Laboratory Q. A de- 
scription of the detectors considered here and of 
the shield can be found in Q . We recall only that 
each detector is viewed by two low-background 
EMI photomultipliers working in coincidence at 
single photoelectron threshold; 2 keV is the soft- 
ware energy threshold. It corresponds in our case 
to > 11 photoelectrons (depending on the detec- 
tor), i.e. well distinguishable from PMT noise. 
In fact, this last is present as fast, single, spare 
photoelectrons, while the "physical" pulses have 
a time distribution with a decay time of hundreds 
ns|. 



10 = 




cl — 














10° 
lO- 2 


















LO 2 


Mill 


1 1 1 1 1 1 
c4 — 


1 1 1 1 1 1 1 1 1 1 1 1 
— c5 — 


1 1 1 1 1 1 1 1 1 1 1 1 










10° 
10-2 










Mill 


1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 


102 
10° 




c7 — 


— c8 — 












10-2 



















5 10 15 20 5 10 15 20 5 10 15 20 
Enerqv CkeV) 



Figure 1. Measured energy distributions for nine 
detectors. 



A statistics of 4549.0 kgxdays is available for 
the annual modulation studies: 3363.8 kgxdays 

1 We comment, in particular, that a pulse shape discrim- 
ination — even under the assumption of an "ideal" elec- 
tromagnetic background rejection — cannot account alone 
for a WIMP signature. In fact, e.g. the neutrons and 
the internal end-range ct's induce signals indistinguishable 
from WIMP induced recoils and cannot be estimated and 
subtracted in any reliable manner at the needed precision. 
2 We recall in addition that scintillators are not affected 
by microphonic noise as ionizing and bolometer detectors. 



in winter time and 1185.2 kgxdays in summer 
period. It is distributed along the year in the 
following way: —1 < cosuj(ti — to) < —0.334 in 
winter time and 0.932 < cosui(ti — to) < 0.996 
in summer time. It is evident that part of the 
time was devoted to systematic studies on cali- 
brations and detector features and qualification 
(more extensively in the winter period, being the 
starting one), that is all the available low en- 
ergy data collected at that time have been con- 
sidered. Obviously the analysis method properly 
takes into account both the collected low energy 
statistics and the cosine range in which it has 
been taken. In fig.l the measured energy distri- 
butions are shown. 

The stability control can profit of the 210 Pb 
peak at 46.5 keV present at level of a few cpd/kg, 
mainly due to a surface contamination by envi- 
ronmental radon during the first period of crys- 
tals storage underground; both the peak position 
and the resolution are controlled binning together 
the data every ~ 7 days (see fig. 2). This allows 
an intrinsic monitoring of the threshold and of 
the calibration stability (e.g. PMT gain and elec- 
tronic line stability). 
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Figure 2. Typical 210 Pb peak at 46.5 keV (see 
text). 



Note that here we do not consider any pulse 
shape analysis (PSA) to reject electromagnetic 
background (see |l]J^]). In fact: i) the PSA has 
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- in any case — a statistical nature; ii) when- 
ever it gives upper limits on the signal (as it is 
generally the case), its results cannot be used for 
the annual modulation analysis; iii) the annual 
modulation analysis acts itself as a very efficient 
background rejection. In any case, as shown in 
fig. 3, the result obtained here is well compatible 
with the upper limits obtained by using PSA con- 
sidering also these detectors and published in [Q . 

3. The maximum likelihood method and 
the quest for a candidate 

To determine the cross section and mass of 
a possible WIMP candidate, a time correlation 
analysis of all the data — properly considering 
the time occurrence and energy of each event — 
has been performed. 

The experimental data collected by a j-th de- 
tector of mass Mj are considered as grouped in 
1-day time bins and in AE = IkeV energy bins; 
the number of the events in the i-th day and k-th 
energy bin, Ny^, will follow a poissonian statis- 
tics with mean value given by fajk = (bjk + 
So,k + S m ,kCosLu(U - to))MjAtiAEejk, where a 
time independent background contribution, bjk, 
has been included in addition to the dark mat- 
ter signal searched for. Here Att represents the 
detector running time during the i-th day (At < 
1 day) and €jk represents the analysis cut effi- 
ciency |lj. The maximum likelihood function re- 
n 3 

suits: L = Tl, ljk e~^^ k . 

Firstly, to allow a direct simple comparison of 
the sensitivity reached here with those of previ- 
ous experiments ^,0, the function y = — 2Zn(L) — 
const = Sfcj/fc, can be minimized with respect to 
(bjfe+So,fc) and S„ l: fc, using for each k-th energy 
bin the data collected at the same time with the 
9 crystals; the y k can be here minimized indepen- 
dently. The obtained S TO) fc values are shown in 
table 1; they point out the relevance of increasing 
the collected statistics. 

However, from this table, a first qualitative 
view can be also obtained. In fact, e.g. in 
the first part of the measured energy region, 2- 
12 keV (where a signal contribution can "a pri- 
on" be expected), < S m > 2 -u= (0.037 ± 0.008) 



Table 1 



S TO values obtained by the maximum likelihood 
method. 



Energy 
(keV) 


(cpd/kg/keV) 


Energy 
(keV) 


(cpd/kg/keV) 


2-3 


0.023 ± 0.037 


11-12 


0.053 ± 0.027 


3-4 


0.017 ± 0.030 


12-13 


0.015 ± 0.028 


4-5 


0.036 ± 0.027 


13-14 


0.017 ± 0.029 


5-6 


0.042 ± 0.021 


14-15 


0.023 ± 0.030 


6-7 


0.038 ± 0.022 


15-16 


-0.045 ± 0.029 


7-8 


0.003 ± 0.023 


16-17 


0.008 ± 0.030 


8-9 


0.050 ± 0.024 


17-18 


-0.031 ± 0.029 


9-10 


0.065 ± 0.025 


18-19 


0.002 ± 0.031 


10-11 


0.032 ± 0.026 


19-20 


0.016 ± 0.030 



cpd/kg/keV [j is obtained, while in the second 
one, 12-20 keV (where any signal is "a priori" ex- 
pected to be suppressed), < S m >l2-20= (0.000± 
0.010) cpd/kg/keV is found. The first value qual- 
itatively supports the presence of a possible mod- 
ulation, requiring therefore a suitable statistical 
analysis of the data to verify if it can be in total 
or in part (as quantified later by the C.L. and the 
X 2 values) ascribed to a possible WIMP presence, 
or it can not. The second value well supports 
the absence of an overall systematics significantly 
exceeding the statistical error in the interesting 
energy region. 

We stress that here and in Q the achieved 
result (see later) is not linked to this qualita- 
tive description (that is, it is not linked to any 
< S rn >2-i2 value), but to the application of the 
maximum likelihood method on all the available 
data (2-20 keV) described in the following. Also 
an unbiased statistical analysis on all the avail- 
able data at the same time is then performed to 
obtain quantitative estimates on a possible effect. 

The maximum likelihood method is well suit- 
able to test the following hypothesis: the modu- 
lated effect qualitatively described above can (or 
not) be accounted by a WIMP with cross sec- 

3 This first qualitative approach is introduced in during 
a preliminary discussion on the Freese et al. method; the 
values were quite similar (see H] for details). 
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tion and mass allowed for istance for the neu- 
tralino |J, assuming a Spin-Independent (SI) in- 
teraction and, if it can, at which C.L.. In partic- 
ular, we refer our results to the quantity £<t p with 
£ = PW p" P , Po = 0.3 GeVcm" 3 and a p WIMP 
cross section on proton. For this purpose, we 
note that So,fc and S mj fc can be written — point- 
ing out their dependence on £a p and the WIMP 
mass, M w — in the form So,k — £,&pSto,k{M w ) 
and S m> k = ^cr p S> mt k{M w ) respectively. The 
S'o,k(M w ) and S/ m ^(M w ) can be calculated ac- 
cording to Q (e.g. besselian form factor, £=1, 
etc.). The £a p and M w values in the best agree- 
ment with the experimental data have been ob- 
tained by minimizing here the y function — using 
all the events with their energy and time occur- 
rence — with respect to the free parameters: £<7 P , 
M m and b^ fc s. A two-step minimization strategy 
allowed us to handle this large number of pa- 
rameters. In fact, by a preliminary y minimiza- 
tion the (bjfc+So,fc)=fjfc values have been deter- 
mined and then, by a subsequent one, also the £a p 
and M^ values. In this last step, the conditions 
M w >25 GeV and (b jfe +S , fc ) = i jk if £op < ^ 
or (bjfc+So,fc) —S.o'pSfQ.k otherwise are required, 
to take into account both the results achieved at 
accelerators for SUSY particles and the obtained 
values for the unmodulated term. 

The minimum value of the y function has 
been found for M w = (59 tT 9 ) GeV and £cr p = 
(1.0 ) 10- 5 pb |. 

4. Consistency checks and statistical eval- 
uations 

To verify the consistency of this result M w has 
been fixed and both £<7 P and the modulation pe- 
riod, T, have been considered as free parameters, 
obtaining still the previous value for £<t p and T 
= (1.3 t° l ) years, a period compatible with a 
yearly modulation. Then, both M w and T have 
been fixed, while £er p and the phase, to, have been 
considered as free parameters, obtaining the al- 
ready found £<7 P value and to = (140 t% ) days, 
a phase compatible with ~ 2 nd June (~ 153-th 
day in the year 1996). The uncertainties on T 
and to are due to the limited available statis- 
tics and to the particular periods of data taking. 



The analysis of the maximum likelihood ratio, 
A = xp^j, has been performed to test the good- 
ness of the null hypothesis Ho (absence of mod- 
ulation) with respect to the Hi hypothesis (pres- 
ence of modulation according to the given M w 
and £<J p ). From the definition of the y function 
we obtain A = e [y(Hi)-y(H )]/2. clearly < A < 1. 

A A value close to 1 will imply absence of modu- 
lation, while a A value close to will support the 
presence of modulation with the given M w and 
£<7 p . To perfom a quantitative statistical test, one 
can use the variable (-2 InA) which is asymptoti- 
cally distributed as a x 2 (high values of -2 InA will 
support presence of modulation). At 90% C.L. 
the upper tail Xo is equal to 2.7 for 1 degree of 
freedom; in our case -21nA = 3.14 > Xq- Therefore 
the hypothesis of no modulation can be rejected 
in favour of the hypothesis of modulation with the 
given M w and £er p at 90% C.L. At this point the 
agreement between Hi and the experimental data 
has been analysed by a x 2 test comparing the ob- 
tained S m of table 1 with the expected values; a 
probability of 6%, to obtain — only because of 
statistical fluctuations — a x 2 value higher than 
we found, has been calculated. This probabil- 
ity is mainly limited by the S m values between 8 
and 12 keV, which show up to ~ 2a fluctuations 
from the expected ones. It is, therefore, evident 
the relevance of data now under analysis with a 
statistics about 7 times larger than the present 
one. We note that here the — widely considered 
— Helm SI form factor has been used for iodine 
|nj; in any case, no relevant effects — within the 
present errors — on the quoted M. w and £<r p were 
observed when adopting different SI form factors 

OS- 

Moreover, we have also analyzed the Nal(Tl) 
data with the procedure described above, to test 
the hypothesis of a modulation effect due to a 
spin-dependent (SD) coupled WIMP. The calcu- 
lation has been performed according to Q, but 
using here for iodine the recently published Ressel 
SD form factor fl2] | . In this way no minimum has 

4 To have a qualitative view one can overimpose the S m 
curve calculated using the found M m and £cr p values on the 
S m of table 1; it is evident that the calculated curve is fully 
under the Sm points and the divergences are quantitatively 
represented by the poor value for the C.L. 



5 



been found for the y function in the considered 
mass region. As mentioned above, preliminary 
results on searches for the annual modulation sig- 
nature have been also obtained by using a liquid 
Xenon scintillator (statistics of 408.2 kgxdays) 
and the Canfranc Nal(Tl) detectors (statistics of 
1342.8 kgxdays) 0. This last experiment had a 
poorer sensitivity than the previous one and, to 
a larger extend, than the present one, mainly be- 
cause of the higher background rate. Also these 
experiments considered only the two extreme pe- 
riods. In particular, to obtain a comparison, we 
have examined the Xenon result of gj. In this 
case we have performed a standard best fit on 
the S m values, constraining M w by the accelera- 
tor limits and £er p by the more stringent results 
obtained with the Nal(Tl) detector Q. An indi- 
cation for M w ~ 60 GeV and £a p ~ 0.5 • 10~ 5 pb 

- values compatible with those obtained above 

— has been found, but with extremely poor C.L.; 
this could be ascribed to the reduced statistics (a 
factor ~ 11 smaller than the present one) and 
sensitivity available there. 




M w (GeV) 

Figure 3. £er p versus M w : the region allowed by 
this preliminary analysis at 90% C.L. is shown 
and superimposed on the best upper limit contour 
for SI interaction, obtained so far [1,13]. 



In fig. 3 the region allowed at 90% C.L. — for 
a SI coupled candidate — by the obtained £cr p 
and M M values is shown; the best upper limit 
contour for SI interaction obtained so far [p],p^[ 
is superimposed. The shaded region represents 
the values for M w and £cr p neither excluded by 
the present analysis nor by the exclusion plot. It 



has been noted that this region is well embedded 
in the Minimal Supersymmetric Standard Model 
(MSSM) estimates for neutralino 

5. Is it a robust case? 

We have always clearly stressed that this pre- 
liminary result needs further investigations [Q, 
however we are frequently addressed with ques- 
tions on its " robustness" . We can comment that 
the present statistical evidence, although in the 
usual range considered in rare event searches, it 
is not very stringent. Only very large exposure 
would possibly allow to reach a firm conclusion; 
similar exposures will be obtained in next future 
by our experiment, which is continuously run- 
ning. 

A long list of "possible" systematics has been 
suggested; however, we have to remark that sev- 
eral of them — if present in an experiment - 
have to be classified as malfunctioning, not as 
systematics. This is the case e.g. of possible un- 
controlled energy threshold and PMT gain vari- 
ations. We stress, in any case, that systematics 
is function of the quality of an experiment, there- 
fore its nature and level is generally very different 
from one experiment to another. On the other 
hand, in the annual modulation case systematics 

— if present — can either simulate the presence 
of an unexisting signal or cancel the presence of a 
real one; therefore, an "a priori" decision on the 
role of its possible " generic" presence is arbitrary. 

In the data taking considered here, we monitor 

— in addition to the controls by energy calibra- 
tion - the external environmental radon, the HP 
N2 flux and the overpressure of the inner Cu box 
(in which the detectors are), the temperature and 
the total and single crystal rates over the single 
photoelectron threshold (i.e. from noise to "in- 
finity"). In fig. 4 examples of the behaviours of 
some of the parameters in long term are shown. 

In the following, we summarize the main fea- 
tures to control systematics: a) data are taken 
first in winter and, then, in summer time (there- 
fore possible positive effects cannot be accounted 
by isotope decay); b) the threshold, the PMT 
gain and electronic line stability are verified both 
by the features of calibration peak and by moni- 
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Figure 4. Example of behaviours of some param- 
eters. From top to bottom: external radon (see 
text); operating temperature; HP N2 flux in the 
Cu box containing the detectors. 



toring the rates; c) the operating temperature is 
controlled and the environmental temperature in 
the installation is not influenced by external sea- 
sonal variations, being conditioned; d) the sta- 
bility in the 12-20 keV energy region allows to 
verify that no appreciable variation of neutron, 
of electromagnetic background and of environ- 
mental conditions were present, although it is not 
clear how all of them could vary with the same 
period and phase of a possible WIMP signal; e) 
the external environmental radon is recorded, al- 
though two levels of sealing in supronyl (main- 
tained in HP N 2 long stored underground) isolate 
the shield containing the Cu box in which the de- 
tectors are closed. The HP N2 flux in the Cu box 
and its overpressure are monitored; f) as regards 
the "possible modulation" spread over the crys- 
tals, for the sake of completeness we recall that 
the S m on single crystals in Q| have been cal- 
culated with another method (see also the com- 
ments there). When properly addressed, obser- 
vations on spread would demonstrate — as first 
- the necessity of long data taking to avoid the 
difficult position of possible "tail" effect (see e.g. 
also other kinds of rare event searches when com- 
paring their results from different periods of data 
taking). In any case, we remark that the C.L. 
found in the overall analysis already accounts for 
the single crystal response (remember the method 



features and the L function definition); g) the ver- 
ification of the period and phase of a possible ef- 
fect by the maximum likelihood method allows 
to restrict mainly to systematics effects able to 
induce an annual modulation with 1 year period 
and ~ 2 nd june phase. Obviously this will be a 
stronger constraint when almost whole years data 
taking (periodical calibrations^] and others are ob- 
viously needed in any case) will be considered. 

Finally, we recall that an almost whole year 
statistics, in improved condition, is already un- 
der analysis and the experiment is continuously 
running, therefore in near future most of uncer- 
tainties will be — in any case — overcome. 

6. Conclusions 

Considering both the difficulty and the rele- 
vance of this kind of searches, a cautious actitude 
is mandatory (see [Q). Here we have presented a 
bare status report on first experimental data and 
analysis method. We are now further carefully 
investigating — with much larger statistics and 
in improved conditions — the region singled out 
in fig. 3. 
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